function F=thbfun(t,u)
global L m k g 
F=[u(2);u(1)*u(4)^2+g*cos(u(3))-k/m*(u(1)-L+m*g/k);
    u(4);-2*u(2)*u(4)/u(1)-g*sin(u(3)/u(1))];
end